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Abstract 

Effects of size and charge asymmetry between oppositely charged ions or particles on spatial 
inhomogeneities are studied for a large range of charge and size ratios. We perform a stability 
analysis of the primitive model (PM) of ionic systems with respect to periodic ordering using the 
collective variables based theory. We extend previous studies [A. Ciach et al., Phys. Rev.E 75, 
051505 (2007)] in several ways. First, we employ a non-local approximation for the reference hard- 
sphere fluid which leads to the Percus-Yevick pair direct correlation functions for the uniform case. 
Second, we use the Weeks-Chandler-Anderson regularization scheme for the Coulomb potential 
inside the hard core. We determine the relevant order parameter connected with the periodic 
ordering and analyze the character of the dominant fluctuations along the A- lines. We show that 
the above-mentioned modifications produce large quantitative and partly qualitative changes in the 
phase diagrams obtained previously. We discuss possible scenarios of the periodic ordering for the 
whole range of size- and charge ratios of the two ionic species, covering electrolytes, ionic liquids, 
charged globular proteins or nanoparticles in aqueous solutions and charge-stabilized colloids. 
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I. INTRODUCTION 



The study of phase diagrams of ionic systems in which the phase separation is mainly 
driven by electrostatic forces is of great fundamental interest and practical importance. Elec- 
trolyte solutions, molten salts, ionic liquids and charge-colloidal suspensions are examples 
of systems with dominant Coulomb interactions. The strong correlations between ions and 
counter-ions are also known to play an important role in determining structure and phase 
behavior of micelles, polyelectrolytes and proteins. Over the last few decades, the phase be- 
havior of ionic fluids has been the subject of many experimental, theoretical and simulation 
studies and the reviews of the state of the art in this field are available in Refs. A 
great deal of the research has been focused on the fluid-fluid phase separation, whereas the 
fluid-solid phase transition has received less attention so far. 

On the other hand, recent experimental, simulation and theoretical studies of room tem- 
perature ionic liquids (RTILs) [7|, charged globular proteins or colloidal particles 
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131 ] reveal structural inhomogeneities of different type and spatial extent, and it be- 



comes evident that the ordering in systems dominated by Coulomb interactions has a very 
rich and complex nature, and is far from being understood. Consequently, ordering of ions 
or charged particles on different length scales attracts increasing attention. 

The much studied simple salts or electrolytes are quite well understood. In these systems, 
the sizes of ions are similar. Therefore, theoretical studies focused mainly on the restricted 
primitive model (RPM) where ions of the same valence are modeled by charged hard spheres 
of the same diameter immersed in a structureless dielectric continuum. At large enough 
volume fractions, transition to an ionic crystal occurs in molten salts when temperature 
decreases 14|-|l7l|. In the ionic crystal the charge is periodically ordered and this ordering 



is subjected to the charge neutrality. On the other hand, charge-ordered, neutral living 
clusters of various sizes and life times were found in simulation studies of the fluid phase 
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2l|. 



A direct inspection of the structure in ionic systems is complicated due to a small size 
of the ions. However, in compliance with the law of corresponding states one my expect a 
similar behavior in oppositely charged colloidal particles of similar sizes under appropriate 
rescaling of the length and energy units. The natural length scale is the sum of radii of 
the anion and the cation, a±, while the energy unit is the Coulomb potential between the 
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oppositely charged ions at contact, Eq = g + g_/(ecr±), where q a is the charge of the a- 
type ion (a = +, — ) and e is the dielectric constant of the solvent. If the non-Coulomb 
interactions are negligible, the phase diagrams of systems consisting of spherical ions of 
various sizes should have the same form in terms of the volume fraction and a reduced 
temperature T* = ksT/EQ. Indeed, simulations show that the phase diagram of charged 
particles with similar sizes in deionized solvent (large Debye screening length) resembles the 



phase diagram of the RPM 
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23jj. In the above reduced units, the room temperature 
is very low for simple salts, whereas for colloidal particles it is high. Ordering of colloids 
at room temperature is directly accessible to observations and vacancies in the crystal can 
be seen [24j. Formation of neutral aggregates in the RPM agrees with recent experimental 
studies of oppositely charged proteins of similar sizes (9), l^j. The aggregation was entirely 
suppressed when the electrostatic interactions had been screened by addition of a sufficient 
amount of salt. Thus, the electrostatic interaction is essential at least for the first step of 
the self-assembly of oppositely charged globular proteins into aggregates. Large spherical 
aggregates presumably indicate nucleation of a crystal. 

RTILs differ from simple salts mainly because the sum of radii of the ions is much larger. 
This reduces the strength of the Coulomb interactions at ion contact and, in turn, increases 
the room temperature in the reduced units. In addition, in many RTILs the size and 
shape of the anion can be quite different from the size and shape of the cation. Moreover, 
non-Coulomb interactions may play an important role here. In such a case the RPM is 
an oversimplification, and the above discussion does not apply to the RTILs. Recent ex- 
perimental and simulation studies on a series of imidazolium-based RTILs indicate spatial 
inhomogeneities on the length scale ~ 10 nm 171. The inhomogeneous structure of liquid 

n 

resembles the structure of the crystal [7j. The origin of such nanoheterogeneity has not been 
convincingly explained yet, and new studies are required. 

In a mixture of positively and negatively charged globular proteins, the size difference 
between the positively and negatively charged molecules can be similar to the size difference 
between the anion and the cation in the RTILs. Recent experimental studies show that 
the size difference plays an important role for the assembly of proteins into aggregates. 
According to Ref. [9|, formation of spherical aggregates of oppositely charged proteins with 
overall charge near zero requires charge- and size compensation. 

The size difference between charged globular proteins and counterions in a solution is 
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much larger than in the previously described cases. For example, the diameter of a lysozyme 
molecule is of order of 3 nm, i.e., it is 10 times larger than the counterions in an aqueous 
solution. The charge of small globular proteins is of order of lOe, where e is the elementary 
charge. Scattering experiments [8:] suggest clustering of protein molecules in water. Effective 
interactions between the lysozyme molecules were shown to have a form of short-range 
attraction and long-range repulsion - the latter resulting from the charges on the molecules 



261 ] . The effective attraction between two molecules is caused by the interactions between 
them and the counterions attracted to both molecules, as well as by the non-Coulomb forces. 
Simulations carried out in Ref. j^j for the potential derived in Ref. j^J show that the fraction 
of clustered molecules and the shape of clusters strongly depend on the lysozyme volume 
fraction. 

Charged nanoparticles or colloidal particles with a diameter ~ 10 4- 1000 nm are several 
orders of magnitude larger than the microscopic counterions in a solution, and their charge 
can be orders of magnitude larger than the charge of the counterion. In these systems, the 
experiments indicate the formation of colloidal crystals with large interparticle distance for 
small volume fractions, and a re-entrant melting for larger volume fractions 1(3, 111. 

From the above description of different experimental systems it follows that spatial in- 
homogeneities are common, but their extent and nature still need to be comprehended and 
classified for different sizes and charges of the anion and the cation. In the first step, it is 
essential to predict the structure formation and phase diagram for various size and charge 
ratios for a generic model where non-Coulomb interactions, deviations from a spherical shape 
and flexibility of molecules are neglected. The role of the above factors can be determined 
after comparing the results obtained within a generic model with the experiments. 

The generic model that allows one to predict the phase separation driven exclusively 
by Coulomb forces is a primitive model (PM). In this model, the ionic fluid is described 
as an electroneutral mixture of charged hard spheres immersed in a structureless dielectric 
continuum. The PM pair potential for two ions a and at distance r apart is 

{oo, r < a a p 
, r > cr a /3 
er 

where an ion of species a has a diameter a a , and charge q a , a a/ 3 = \ [a a + ap) and e is the 
dielectric constant. The PM is the simplest model for all the above discussed systems. The 



4 



two-component PM can be characterized by the parameters of size and charge asymmetry: 

A = ^±, Z=^. (2) 
°- ' \Q-\ 

For X = Z = 1 one arrives at the RPM. 

In the PM with large difference in the ion radii, the spatial distribution of the ions is 
expected to be quite different from the RPM, because the tendency for minimizing the 
electrostatic e nerg y competes with the geometrical restriction on packing of spheres with 
different sizes 24|. Packing of large and small spheres that maximizes the entropy could 
lead to the formation of mesoscopic charged regions, and, on the other hand, when the 
electrostatic energy is minimized, periodic pattern involving voids could be formed. As 
a result of the competition between maximizing entropy and minimizing energy, both the 
charge and the number- density of ions can oscillate in space. 

The systematic studies of the effects of size and charge asymmetry on the periodic ordering 
in PMs were initiated in Ref. j^sj within the framework of the field theoretical description. 
Based on a mean-field stability analysis, the authors found the boundaries of stability of the 
disordered phase for the whole range of A and Z. It was shown that, besides a gas-liquid 
separation, in a certain portion of the phase diagram, the uniform fluid became unstable 
with respect to the order parameter oscillations of wavelengths 2ir/kf) with kb ^ 0. The line 
in the phase diagram corresponding to the instability of the disordered phase with respect 



to periodic ordering is called the A-line 
kb = 0. The results obtained in Ref. 



29| to distinguish it from the spinodal line for which 



show that (i) the periodic ordering mainly depends 
on the size asymmetry; (ii) the qualitative dependence on the charge asymmetry is found 
only for a sufficiently large size asymmetry. 

In this paper we continue the systematic study of the periodic ordering in asymmetric 
PMs. We extend the previous study in several ways. The first modification concerns an 
approximate description of the reference hard-sphere mixture. In Ref. [28j a local-density 
approximation is employed for a hard-sphere free energy functional. Here, we consider 
a non-local approximation for the reference hard-sphere fluid which is shown to lead to 



30] . In our calculations we use the 



the Percus-Yevick (PY) theory for the uniform case 
Lebowitz's solution of the generalized PY equation [31]. 

Another modification concerns the regularization of the Coulomb potential inside the hard 
core. It is worth noting that in the treatments of models with hard cores, the perturbation 
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potential is not denned uniquely inside the hard core. Here we use the Weeks- Chandler- 



Andersen (WCA) regularization scheme for the Coulomb potentials 0^/9 ( r ) 



er 



r < a al3 



(3) 



instead of 0?g(V) = q a qp8( r — Cap)/( er ) adopted in 28]. As was shown in Ref. |33|, the 
simple form for </>^g(r) given in Eq. ([3]) produces rapid convergence of the series of the 
perturbation theory for the free energy. On the other hand, the best theoretical estimates 
for the gas-liquid critical point of the RPM was obtained using the WCA regularization 



scheme 



34]. 



Finally, we extend the study of the relevant order parameter (OP) undertaken in 



Following the ideas of Refs. [35|,|36], we determine the OP connected with the phase transition 
to an ordered phase and analyze the character of the dominant fluctuations along the A-lines 
associated with the periodic ordering. 

A theoretical background for this study is the statistical field theory that exploits the 
method of collective variables (CVs) 37H40|. The theory enables us to derive an exact 



expression for the functional of grand partition function (GPF) of the model and on this 
basis to develop the perturbation theory 4CN42l|. As was shown in Ref. 42j, the well-known 
approximations for the free energy, in particular Debye-Huckel limiting law and the mean 
spherical approximation, can be reproduced within the framework of this theory. Links 
between this approach and the field theoretical approach [28] were established in [43j for the 
case of the RPM. 

Our paper is organized as follows. In Section 2 we give some brief background to the 
CVs based theory for the PM. Based on the Gaussian approximation of the functional of 
GPF we obtain the pair direct correlation functions and determine the OP characterizing 
the periodic ordering in PMs. In Section 3 we study the effects of size and charge asymmetry 
on the periodic ordering taking into account the above-listed modifications. We discuss the 
results in Section 4 and conclude in Section 5. 
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II. THEORETICAL BACKGROUND 



A. Functional representation 

We start with the general case of a two-component PM consisting of N + cations carrying 
a charge q + = Zq of diameter a + and iV_ anions carrying a charge g_ = — q of diameter o~_. 
The ions are immersed in a structureless dielectric continuum. The system is electrically 
neutral: J2 a =+ - QaPa — and p a = N a /V is the number density of the ath species. 

The pair interaction potential is assumed to be of the following form: 

U a p{r) = <Q{r) + 4>%{r), (4) 

where </>?f (r) is the interaction potential between the two additive hard spheres of diameters 
cr a and a p. We call the two-component hard-sphere system a reference system. Thermody- 
namic and structural properties of the reference system are assumed to be known. 09g(r) is 
the Coulomb potential given in Eq. d3J), and hereafter we put e = 1. 

Using the CV method we get an exact functional representation of GPF for the PM with 
size and charge asymmetry |4o| 



I (dp)(dw)exp ( -— ^2^2^(k)p k>a p^+ 

\ a,/3 k 

+1 ^2 W k,QPk,a + hi E m [u a - iu a ] J . 
a k / 



(5) 



In Eq. (jSJ) Pk,a = Pk,a — iPk,a i s the CV which describes the value of the k-th fluctuation 
mode of the number density of the ath species, each of p k a (p k a ) takes all the real values 
from — oo to +oo; w kja is conjugate to the CV p kjCl ; (dp) and (da;) are volume elements of 
the CV phase space 

( d P) = II dPo,«n dPk, Q dp kia , (dw) = JJ dwo.aJJ dw ka dc4 Q 

a k^O as k^O 

and the product over k is performed in the upper semi-space (p-k.a = Pk a , w -k,a = wL)' 
4>ap(k) is the Fourier transform of the Coulomb potential. In the case of the WCA 
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regularization (see Eq. ([3])) we obtain for P^Jk) |41| : 



T*(l + 5) x 



T*Z{\-5) x 
««_<*) = -*g*E&, (8, 



where the following notations are introduced: 



is the dimensionless temperature, x = fccr±, o± = (cr + + cr_)/2 and 

„ A-l 



A+l 

Similarly, hereafter we introduce the parameter v 

Z- 1 



(10) 



(11) 



Z + l 

The parameters 5 and v are more convenient than the parameters A and Z because they 



vary between —1 and 1. Following Ref. 



28|, we choose the dimensionless temperature T* 



given in ([9]) and the volume fraction of all ions 



C = ^(p + 4 + P-^ 3 ) (12) 



as thermodynamic variables. 

Shs[^o — i^a] is the GPF of a two-component hard-sphere system with the renormalized 
chemical potential 

k 

in the presence of the local field —iu a (r). In ( !T3|) u a is the dimensionless chemical potential, 
v a = (3fi a — 31nA a , \x a is the chemical potential of the ath species, f3 is the reciprocal 
temperature, A" 1 = (2Trm a /3^ 1 / h 2 ) 1 ^ 2 is the inverse de Broglie thermal wavelength. 

In order to develop the perturbation theory we present In Shs[^q — iw a ] in the form of the 
cumulant expansion 

-ir 



IhShs[...] = X] — J yX ai ... an (v a ;ki,...,k n ) x 

n>0 ai,...,a n ki,...,k n 

xu ki,«i • • • w kniQn 5 kl+ ... +kn , (14) 



where 5ki+...+k n is the Kronecker symbol. In Eq. ([H]) the nth cumulant Wt ai ... an coin- 
cides with the Fourier transform of the n-particle connected correlation function of a two- 



component hard-sphere system 40]. 



B. Gaussian approximation 

Having set 0Jt Q , 1 ... Q , n = for n > 3, after integration in Eq. §5§ over oj^ a one arrives at 
the Gaussian approximation for the functional of GPF 41] 

E G [u a ] = E MF [u a ] 5' / (dp)exp | - — ^^C a ^(A;)/) kiC( p_ k ^|, (15) 

a,/3 k 

where E MF is the GPF in the mean-field approximation and S' = Y\ det[\^ 2 9Jt 2 ]~ 1//2 with 

k 

VJt 2 being the matrix of elements 9Jt Q( g(/c) /V '. C a p(k) is the Fourier transform of the pair 
direct (vertex) correlation function in the random phase approximation (RPA) 

C a0 (k) = P4>^(k)+C™(k). (16) 

In ( I16p . C^(k) is the Fourier transform of the pair direct correlation function of a two- 
component hard-sphere system. It is connected with ^K a p{k) by the relation C^ s (k)DJl2{k) = 
1, where C^ s (k) denotes the matrix of elements C^§(k) and 1 is the unit matrix. In the 
limit k = 0, the CSf coincides with the coefficients a a p which are obtained in j^sj as a result 
of the local-density approximation. 

It is convenient to introduce CVs which describe the fluctuation modes of the total number 
and charge density, pk,7v and Pk,Q, by the relations 

pkn = y+z (pk,+ + pk,-) , 

Pk,Q = z^ r ^(Zp k ,+ - fr-). (17) 
Then, Eq. (ITS"]) can be rewritten in terms of p k ,Af an d Pk,Q as follows: 

E G [v a ) = EuF^a] 2' J (dpAr)(dp Q )exp| - ^-^2^2c A B(k)p KA p^ KB y (18) 

A,B k 
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where A(B) = N,Q and 

C QQ (k) 
C QN (k) 



1 

TTzy 
1 



C++(k) + Z 2 C__(k) + 2ZC + _(k) 
C++(k)+C—(k)-2C+-(k) 
C ++ (k) - ZC^(k) + (Z- l)C+-(k) 



(19) 



(i + zy 

are the density- density, charge-charge and charge- density direct correlation functions, re- 
spectively 

In general, an equation for the boundary of stability of the uniform phase with respect 
to fluctuations is given by 



det C 5 



k=kh 



o 5 



(20) 



where Ci denotes the matrix of elements Cab(&) (or C a p(k)). The corresponding wave vector 
kb is determined from the equation [29] 

<9det C 2 /<% = 0. (21) 

The case k = kb = corresponds to the gas-liquid-like separation 28|, |41|. Here we are 
interested in the A-line, the boundary of stability associated with fluctuations of the OP 
with k = kb 7^ 0. On the A-line, the fluid becomes unstable with respect to the periodic 
ordering indicating that there can be a phase transition to an ordered phase. 



C. Order parameter 



The determination of the OP is the important issue in the phase transition theory of 
mixtures. This problem has got a consistent and clear solution within the given approach. 
In order to determine the OP associated with the periodic ordering we follow the ideas of 



Refs. 



35 



361 ] . First, we diagonalize the square form in Eq. (118p by means of an orthogonal 
transformation 

£k,l = tNNP]s.,N + tNQPk,Q, (22) 

£k,2 = tQNPk,N + tQQpk,Q- (23) 

The explicit expression for coefficients t^B are given in Appendix. The corresponding eigen- 
values Ei(k) and E2(k) are found to be 



ei, 2 (A0 = 2 ( C NN(k) + C QQ (k) ± (C NN (k) - C QQ (k)Y + AC 2 NQ (k) 



1/2 



(24) 
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Based on the solutions of equations (l20l) - (!2l!) we find the eigenvalue which becomes equal 
to zero along the calculated A-lines for the fixed values of parameters 5 and v. We suggest 
that the corresponding eigenmode is connected with the relevant OP. 




FIG. 1: PM with a small size asymmetry (5 = 0.2): the dependence of eigenvalues £\ and £2 on 
the wave numbers (x = ka±). Solid and dashed lines correspond to v = 0.33 (£ = 0.1, T* ~ 0.033) 
and v = 0.9 (C = 0.1, T* ~ 0.068), respectively. 



Our analysis shows that only one of the eigenvalues, i.e., Ez(h), becomes zero along the 
Adine of the both symmetric and asymmetric PMs. The same result was obtained earlier 



for a mixture of neutral particles 



35 



361 ] . In Figs. 1 and 2 we show a typical behaviour of 



£i(k) and £2(^0 under thermodynamic conditions corresponding to the points located on the 
Adines associated with k b 7^ 0. 

We assume that CV £k,2 is connected to the relevant OP. Based on Eqs. (12"2]) -( 123]) we can 
determine the direction of strong fluctuations along the A-line by the relation 



tan# 



(25) 



Inn tQQ 

where 9 is the rotation angle of axes £k 6 ,i and £1^,2 in the plane (pk b ,N,Pk b ,Q)- The case 
9 = corresponds to the pure charge density fluctuations (axes £k 6 ,2 an d Pk 6 ,Q coincide) 
and the case 9 = =F7r/2 corresponds to the pure total number density fluctuations (axis £k 6 ,2 
coincides with axes ±pk,,,7V; respectively). Taking into account the formulas from Appendix 
(Eqs. ((27D-03OD), we have 

tan# = — - = a 1 . (26) 

«2 
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It is worth noting that in the long-wavelength limit a 2 (k = 0) = 0, and one gets 9 = —tt/2 in 
agreement with the expected separation into homogeneous charge neutral dilute and dense 
phases. 
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FIG. 2: PM with a moderate size asymmetry (5 = 0.5, v = 0.33): the dependence of eigenvalues 
E\ and £2 on the wave numbers {x = ka±). Solid and the dashed lines correspond to ( = 0.1, 
T* ~ 0.033 and ( = 0.3, T* ~ 0.037, respectively. 



In Ref. [28|| the analysis was similar, except that the eigenvalues denoted by C^ik) and 
Cr]r](k) were defined directly in terms of C a p in such a way that e 2 (k) = C^k) when 
C+-(k) > and e 2 (k) = C m (k) when C + -(k) < 0. C^{k) and defined in Ref. [28] 

reduce to the direct correlation functions for the charge and for the number density, respec- 
tively for Z = 1 (i.e. for the RPM), and the corresponding eigenmodes reduce for Z = 1 to 
the charge and the number density waves. The advantage of the present approach is that 
the critical mode is associated with the same eigenvalue £2{k) independently of k, ( and T. 
It is also more convenient to present the nature of the eigenmode in terms of the angle 8, 
rather than in terms of the parameter R introduced in Ref. 28( in a way analogous to tan 9 
in Eq.([25]). 



III. BOUNDARY OF STABILITY ASSOCIATED WITH PERIODIC ORDERING: 
RANDOM PHASE APPROXIMATION 

In this section we study the effects of size and charge asymmetry on the boundary of 
stability against the fluctuations with k = k b ^ 0. To this end, we use Eqs. (T20]) - (|2X]) taking 
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into account Eqs. ©-([8]). In order to determine the character of the dominant fluctuations 
we use Eq. fl25|) . 

We take into account the ^-dependence of the direct correlation functions of the refer- 
ence system using an exact solution of the generalized PY equation obtained by Lebowitz 
3jJ . The explicit expressions for the Fourier transforms of the OZ partial direct correlation 



functions of a two-component hard sphere system are given in Refs. 
cumbersome to be reproduced here. 



44 



451 ] . They are too 



As in Ref. 



281 ] . we distinguish the three regimes in size asymmetry: small size asymmetry, 



moderate and large size asymmetry, and very large size and charge asymmetry. Each regime 
is characterized by a typical behaviour of the A-lines. As we will see below, the 5-ranges 



of these regimes slightly differ when compared to 
symmetric case. 



281 ] . We also dwell briefly on the size- 



A. Size-symmetric PM 



We start with a size-symmetric PM corresponding to 5 = (or A = 1). In Fig. 3 the A-line 
associated with the wave vector kh ^ is displayed. As is seen, it is a straight line identical 



to that obtained for the RPM 
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47j: T*(k = k b ) = S x (. For the WCA regularization, 



S x ^ 0.285 and x b = k b a± = k b a ~ 4.078 [47]. 




FIG. 3: The A-line for the transition to the ordered phase for 5 = 0. Temperature T* and the 
volume fraction of ions £ are in dimensional reduced units defined in Eqs. ([9]) and (|12p. respectively. 



For the size-symmetric PM, the angle 9 indicating the direction of the strong fluctuations 
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is equal to zero. Therefore, the A-line shown in Fig. 3 is the boundary of stability of the 
uniform phase against the charge density fluctuations. 
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FIG. 4: The A-lines for the transition to the ordered phase for 5 = 0.1 (y > and v < 0). Lines 
from the top to the bottom: v = 0.9, v = 0.67, v = 0.33, v = 0, v = —0.33, v = —0.67, v = —0.9. 
Temperature T* and the volume fraction of ions C are i n dimensional reduced units defined in 
Eqs. ([9]) and (fl2j) . respectively. 



The charge asymmetry has no effect on the A-line of the size-symmetric PM. This property 
does not persist if the higher-order terms are taken into account in Eq. (fl~4l) (see Refs. [42[ 
Q, S])- The effects of char ge-d ensity fluctuations on a phase behaviour of RPM (5 = 0, 
v = 0) were studied in Ref. [15| using the field theoretical description. It was found that 
in the presence of fluctuations, the A-line disappears. Instead, a fluctuation-induced first- 
order transition to an ionic crystal appears. Interestingly, the wavelength of the charge 
wave is independent of the density along the liquid- crystal coexistence line [151 ] . which was 
interpreted as the formation of vacancies in the crystal of a fixed unit cell when the density 
of ions decreases. We expect the similar situation to take place for the PM with the small 
size asymmetry. 



B. Small size asymmetry 



Now we consider the asymmetric PMs with a small size asymmetry (5 < 0.3 or A < 2). 
This case corresponds to molten salts, electrolytes, some RTIL and to oppositely charged 
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globular protein or nanoparticle mixture in the limit of infinite screening length. The A-lines 
associated with the wave vectors kj, 7^ are shown in Fig. 4 for 5 = 0.1. It should be noted 
that the A-lines are located at the temperatures that are by an order of magnitude lower 
than those in Ref. j^. In addition, the monotonously increasing behaviour of the A-line 
temperature, T* versus (, is found for a sufficiently large charge asymmetry > 0.33). 
For \v\ < 0.33, T* has got a maximum in the range ( ~ 0.42 -j- 0,48. The effect comes into 
prominence with an increase of 5. This differs from the previous results [28] demonstrating 
a monotonous increase of T* with ( along the A-lines (( = 0-r0.7) for all v. The difference is 
directly related to the non-local approximation for the reference hard-sphere system adopted 
in the present work. As is seen from Fig. 4, the A-line temperature increases when v increases 
from —0.9 to 0.9 which qualitatively agrees with 



4.8 



4.4 



4.0 



3.6 



3.2 



0.0 




0.1 



0.2 0.3 



0.4 



0.5 



FIG. 5: The wave number = h^a± corresponding to the ordering of ions along the A-lines shown 
in Fig. 4 for 6 = 0.1 [y > and v < 0). Lines from the bottom to the top: v = 0.9, v = 0.67, 
v = 0.33, v = 0, v = —0.33, v = —0.67, v = —0.9. C is the volume fraction of ions. 



For a small size asymmetry, the wave numbers characterizing the period of the OP oscil- 
lations in a non-uniform phase are x\, = kb<J± > it and their magnitudes depend very slightly 
on ( (see Fig. 5). The comparison with Ref. {28 1 implies that the magnitude of Xb is mainly 
determined by the regularization method of the Coulomb potential inside the hard core. In 
Fig. 5 we also demonstrate the effect of charge asymmetry on x^. As is seen, x^ increases 
with the variation of v from 0.9 to —0.9 for the fixed ( which agrees with the results obtained 
in \2m. 
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FIG. 6: The variation of the angle 9 along the A- lines presented in Fig. 4 (y > and v < 0). 6 is 
measured in degrees. 



We have calculated the angle 9 showing the direction of strong fluctuations along the 
A-line. The results for 5 = 0.1 are presented in Fig. 6. As is seen, 9 changes around zero and 
this change increases with an increase of size asymmetry. As regards the charge asymmetry, 
the modulus of 9 decreases with the variation of the charge asymmetry parameter v from 0.9 
to —0.9 for the fixed S, and \9\ approaches zero for v = —0.9. In particular, for 5 = 0.1 the 
angle 9 changes continuously in the range from —14° to +16.5° for v = 0.9 and from —0.3° 
to +0.75° for v = -0.9 when ( is varied from to 0.5. For S = 0.2 we have -30° < 9 < +44° 
{y = 0.9) and —0.03° < 9 < +1.3° (y = —0.9) in the same range of (. In general, the charge 
density fluctuations are the dominant fluctuations for 5 < 0.3. It should be noted that for a 
small size asymmetry, the A-lines associated with the periodic ordering are located at much 
higher temperatures than the spinodals indicating phase separation into two uniform phases. 

Summarizing, in the PM with a small size asymmetry we expect the phase transition to 
an ionic crystal with a compact unit cell where nearest neighbours are oppositely charged, at 
least when ( is sufficiently large. In particular, the CsCl crystal was observed experimentally 
for the system of oppositely charged colloids of comparable sizes 24| . For low volume fraction 
compact charge-ordered clusters are expected, in agreement with experiments for oppositely 
charged proteins |9| and with simulations 25]. 
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C. Moderate and large size asymmetry 



Let us consider the case of moderate and large size asymmetry corresponding to 5 > 0.3 
(A > 2). In fact, PMs with S = 0.3 demonstrate a crossover-type behaviour. The A-lines 
and the corresponding wave numbers for 5 = 0.3 are shown in Figs. 7 and 8, respectively. 
It follows from our calculations that PMs with S = 0.3 and v < —0.67 do not undergo an 
instability with respect to fluctuations with fc& ^ 0. This is contrary to the results obtained 
for PMs with a small size asymmetry (see Fig. 4). Similar to the small size asymmetry case, 
the A-line temperature decreases when v varies from 0.9 to —0.67. For v = —0.67, the A-line 
associated with ^ is located at T* slightly lower than T* at the spinodal indicating the 
separation in two uniform phases for the same (. 




0.0 0.1 0.2 0.3 0.4 



FIG. 7: The A-lines for the transition to the ordered phase for 5 = 0.3 (y > and v < 0). 
Temperature T* and the volume fraction of ions £ are in dimensional reduced units defined in 
Eqs. ([9]) and (fT2]) . respectively. The curve for v = —0.67 is located at a very low temperature and 
is indistinguishable in this plot. 

As is seen from Fig. 8, the dependence of the wave numbers Xb on the volume fraction £ 
is more prominent than that shown in Fig. 5. In particular, Xf, is an increasing function of ( 
for \u\ < 0.33. Such behavior is consistent with the decrease of the interparticle separation 
for decreasing average volume per particle for the same type of structure. For v > 0.67, Xf, 
first very slowly increases and then again slowly decreases when volume fraction increases. 
For v = —0.67, Xb rapidly decreases. For very small values of (, Xb ~ 7r (lb ~ er+ + a J) 
without regard to the charge asymmetry. + 
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0.1 0.2 0.3 0.4 0.5 



FIG. 8: The wave number Xb = kba± corresponding to the ordering of ions along the A-lines shown 
in Fig. 7 for 5 = 0.3. £ is the volume fraction of ions. 

The variation of the direction of strong fluctuations along the A-lines (angle 9) is shown in 
Fig. 9. As is seen, the charge density fluctuations still prevail over the total number density 
fluctuations for —0.33 < v < 0.9 except for a high-density region (see the case v = 0.9). 

The fluid-solid phase coexistence in the PM with 5 = 0.3 and v = was studied by 



computer simulations in Ref. 50J. In particular, it was found that there is a coexistence at 
T* = 0.033 between a fluid at £ = 0.38 and a (NaCl) solid phase at 0.51. Our calculations 
of the A-line temperature corresponding to ( = 0.38 yield T* = 0.037 (see Fig. 7). At 
this point, the order parameter is connected to the charge density fluctuations (6 ~ 0°) 
suggesting the phase transition to the ionic crystal. 

For 5 > 0.3 all the A-lines associated with kb ^ demonstrate a single maximum for 
the volume fraction ( m (see Figs. 10, 11 and 14). This implies that the periodic orde ring 



is less favorable at higher volume fractions. This behaviour is similar to that found in 28] 



but for the higher values of 5, namely for S > 0.4. It is worth noting that the A-lines lie 



at temperatures lower than those found in 28|. This tendency is kept for the whole region 



of the variation of 5. On the other hand, the comparative analysis performed within the 
same regularization scheme, i.e., the WCA regularization, shows that the account of the 
/c-dependence of the reference system correlation functions leads to a pronounced increase of 
the A-lines temperatures for a moderate and large size asymmetry (see Fig. 10). Such strong 
impact of the above-mentioned fc-dependence on the magnitude of the A-lines temperature 
is contrary to the results found for a one-component fluid with competing attractive and 
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repulsive interactions 
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FIG. 9: The variation of angle 8 along the A-lines presented in Fig. 7. 9 is measured in degrees. 
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FIG. 10: The A-lines for the transition to the ordered phase for v = 0.9 and for two values of 
the size asymmetry (<5 = 0.5 and 5 = 0.6). Solid and dashed lines correspond to a non-local and 
local approximations for the hard-sphere reference system, respectively. Temperature T* and the 
volume fraction of ions £ are in dimensional reduced units defined in Eqs. ([9]) and (|12j) . respectively. 



Unlike Ref. 28], our results demonstrate the dependence of ( m on the both parameters 5 
and v. In general, ( m decreases noticeably with an increase of 5 and has a non-monotonous 
behavior with the variation of v. By contrast, the maximal value of the A-lines temperature 
slightly depends on the size asymmetry for 5 > 0.3. 

Let us describe in some detail the case 5 = 0.5 (A = 3, volume ratio ~ 30), characterizing 
moderate size asymmetry. 5 = 0.5 can be found in a solution of charged molecules of a 
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diameter ~ 1 nm in a presence of counterions of a diameter ~ 0.3 nm, or in some RTIL and 
finally in a mixture of oppositely charged globular proteins or nanoparticles in deionized 
solvent. 
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*H °< 15 
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0,0 0,1 0,2 0,3 0,4 0,5 

FIG. 11: The A-lines for the transition to the ordered phase for 5 = 0.5. Temperature T* and the 
volume fraction of ions £ are in dimensional reduced units defined in Eqs. ([9]) and (|12p . respectively. 

For 5 = 0.5, the wave number Xb increases along the A-line from Xb < vr for Q < 0.2 to 
Xb > vr for ( > 0.2 except for the case of large charge asymmetry v = 0.9 (see Fig. 12). 
Such a behaviour implies that different ordered structures can be formed with the variation 
of (. For v = 0.9, the wave number Xb < vr along the A-line (0 < ( < 0.5) and depends 
only slightly on the volume fraction. For |z/| < 0.33, Xb change their trend sharply for the 
large values of ( (( ~ 0.4). The dominant field is pk,/v at small ( for all v (see Fig. 15). 
For v > 0.33, the character of the dominant fluctuations changes with an increase of (. In 
particular, 9 varies from —61° to +71° for v = 0.9 and from —62° to +49° for v = 0.8 
when the volume fraction increases from to 0.5. For v = 0.67, the dominant field is pk,<? 
when ( > 0.11. For \u\ < 0.33, the angle 9 takes a maximum value 9 m for ( ~ 0.4 and the 
modulus of 9 m increases when v varies from 0.33 to —0.33: mod 9 m ~ 22° for v = 0.33 
and mod 9 m ~ 38° for v = —0.33. Thus, the charge density fluctuations prevail over 
the total number density fluctuations along the Adines for this range of v. In general, for 
the intermediate size asymmetry both the charge density and the total number density are 
inhomogeneous in space. 

Let us focus on the large size asymmetry, 5 > 0.6 (A > 4, volume ratio > 60). For 
5 > 0.6, we obtain Xb < vr in the region < ( < 0.5. The Adine and the dependence of Xb on 
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FIG. 12: The wave number Xb = ki,a± corresponding to the ordering of ions along the A-lines 
shown in Fig. 11 for 5 = 0.5. £ is the volume fraction of ions. 
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FIG. 13: The variation of angle 6 along the A-lines presented in Fig. 11. 9 is measured in degrees. 

the volume fraction for 5 = 0.8 are shown in Figs. 14 and 15. As is seen, Xb first increases 
reaching its maximal value and then it decreases with the increase of (. The maximum 
becomes more pronounced with the decrease of charge asymmetry. Such a behaviour of Xb 
is typical of the systems with 0.6 < 5 < 0.9. The variation of the angle 9 along the A-lines 
for 5 = 0.8 is shown in Fig. 16. For the charge asymmetry range < v < 0.8, the total 
number density fluctuations are the dominant fluctuations along the A-lines. Moreover, the 
contribution from this type of fluctuations increases when the charge asymmetry decreases. 
For v = 0.9, the character of the dominant fluctuations changes continuously along the 
A-line from the total number density fluctuations for < ( < 0.12 to the charge density 
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FIG. 14: The A-lines for the transition to the ordered phase for 5 = 0.8. Temperature T* and the 
volume fraction of ions £ are in dimensional reduced units defined in Eqs. ([9]) and (|12p. respectively. 




FIG. 15: The wave number Xf, = kt,a± corresponding to the ordering of ions along the A-lines 
shown in Fig. 11 for 5 = 0.8. Q is the volume fraction of ions. 

fluctuations reaching a maximal value at £ ~ 0.2. For ( > 0.28, the dominant fluctuations 
are again the total density fluctuations. It should be noted that for 5 = 0.9 the total density 
fluctuations play the dominant role along the whole A-line for < v < 0.9. 

In the approximation considered, PMs with 5 > 0.6 and v < do not undergo instabilities 
against the fluctuations with fc& 7^ 0. Moreover, the A-lines tend to lower T* with the decrease 
of the charge asymmetry coming close to the spinodals associated with the phase separation 
into uniform phases. This implies that the phase transition to the periodic ordering in 
the PMs with large size asymmetry becomes less favorable when the charge asymmetry 
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decreases. 

In conclusion, for the PMs with moderate and large size and charge asymmetry, we expect 
the phase transition to the colloid crystals of different structure formed by a nearly charge- 
neutral units for small and large volume fractions. The size of these units depends on the 
charge ratio and is generally larger than in the case of 5 < 0.3. For volume fractions near 
the maximum temperature at the A-line the number density in positively charged regions 
is different from the number density in negatively charged regions in the fluctuation that 
destabilizes the homogeneous phase (see Figs. 14 and 16). Global stability, however, may 
correspond to another structure. For a large size asymmetry case, the complex structures 
discussed above may be preempted by the phase separation in two uniform phases when the 
charge asymmetry decreases. 
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FIG. 16: The variation of angle 6 along the A-lines presented in Fig. 11. 9 is measured in degrees. 
D. Very large asymmetry 

Let us consider the system with a very large asymmetry in size and charge: 5,v> — > 1. 
For 5 = 0.99 (A = 200) and v = 0.99 (Z = 200), the A-line and the corresponding wave 
vectors are shown in Figs. 17 and 18. As expected, the region bounded by the A-line extends 
to higher temperatures in this case. The A-line assumes a maximum with the coordinates 
( m ~ 0.09, ~ 1.4 and kb, m &± — 2. Remarkably, is five times higher and ( m is two 
times lover than for 5 — v — 0.9 (A = Z = 19). On the other hand, the above-mentioned 
value of is by an order of magnitude smaller than that obtained for the same model in 
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Ref. [28|| . Recall that ( m has no size and charge dependence in the approximation used in 
the previous work. The period of the OP oscillations lb = 27i/ki, ~ na + /2 persists in the 
region < ( < 0.25 and increases for a higher value of the volume fraction. The behaviour 
of the angle 9 indicating the direction of the dominant fluctuations is shown in Fig. 19. As 
is seen, 9 changes its trend sharply from ~ —70° to ~ 70° in the region < C, < 0.05. Then, 
9 remains nearly constant (9 m 80°) in a wide region of ( (0.05 < ( < 0.25) and again 
changes sharply to ~ —90° keeping this value for ( < 0.5. We suggest that such a behaviour 



of the O 

systems 



indicates the re-entrant phase transition observed experimentally in the colloidal 
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FIG. 17: The A-lines for the transition to the ordered phase for 6 = v = 0.99. Temperature T* 
and the volume fraction of ions £ are in dimensional reduced units defined in Eqs. ([9]) and (|12p . 
respectively. 



IV. DISCUSSION 

The present analysis has shown that the size and charge asymmetry is expected to have 
a significant effect on the periodic ordering in the systems with dominant Coulomb interac- 
tions. 

For the PM with a small size asymmetry (A < 2), the trend of the A-lines in a wide region 
of the volume fractions is qualitatively similar to that found for the RPM . This is consistent 
with the behavior of OP demonstrating only a small deviation from the pure charge density 
oscillations. In this case, we expect that the ordered phase corresponding to an ionic crystal 
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FIG. 18: The wave number Xb = kba± corresponding to the ordering of ions along the A-line shown 
in Fig. 17. £ is the volume fraction of ions. 
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FIG. 19: The variation of the angle 9 along the A-lines shown in Fig. 17. 9 is measured in degrees. 



with a compact unit cell is formed at lower temperatures and at higher volume fractions 
when compared to the A-lines. In particular, such a behaviour was confirmed experimentally 
for the system of oppositely charged colloidal particles of comparable sizes (see Ref. 24j). 
On the other hand, we expect that charge-ordered living clusters are formed in the fluid 
phase with the period of the charge wave similar to the one in the crystal. It is worth noting 
that in the case of moderate and small charge asymmetry, the trends of the A-lines deviate 
from the RPM-like behaviour when the volume fraction increases. 

PMs with A = 2 show a crossover-type behaviour between the regime of small size asym- 
metry and the regime of moderate and large size asymmetry. In qualitative agreement with 
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the results of Ref. 



50|, our results predict the phase transition to the ionic crystal for the 



system with A = 2 and Z — 1. 

For a moderate and large size asymmetry (2 < A < 20), the A- lines show a single max- 
imum at volume fraction ( m . While depends slightly on the size asymmetry, ( m 
decreases noticeably when the size asymmetry increases. In addition, decreases with 
the decrease of charge asymmetry. The nonmonotonic temperature at the A-line indiates 
that some kind of a crystalline phases can be stable for intermediate volume fractions, and 
re-entrant melting occurs at the high volume fractions (where the temperature at the A-line 
decreses). It should be noted that both the charge density and the total number density 
oscillate along the A-lines and their contributions to the OP vary depending on the charge 
asymmetry as well as on the volume fraction. Note that for overall charge neutrality the 
number of negatively charged species is Z times larger than the number of positively charged 
species. For a large charge asymmetry the main contribution to the number density comes 
from the negatively charged species. 

The expected phase transition to the crystal phase of different structure is in qualitative 



of large and small oppositely charged spherical colloids with A ~ 3 
three of the predicted structures were also observed experimentally 
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24 1 . Remarkably, 



24j. Our results 



also show that the large nearly neutral clusters may be formed in the PMs with a moderate 
and large size asymmetry at the small volume fractions. Such a behavior was observed 
experimentally in water solutions of the charged globular proteins |Sj. On the other hand, the 
appearance of instabilities in the fluid phase against the periodic ordering for a moderate size 
and charge asymmetry may explain to some extent the inhomogeneous structure observed 
in RTILs. In these cases the role of non-Coulomb interactions should be clarified in future 
studies. The trend of the A-lines found for the large size asymmetry implies that the periodic 
ordering becomes less favorable when the charge asymmetry decreases and the fluid-solid 
phase transition may be preempted by the phase separation in two uniform phases when the 
charge asymmetry decreases. 

For a very large size and charge asymmetry (A = Z — 200) the number density fluctu- 
ations dominate except for small ranges of ( (Fig. 19). The period of the OP oscillations 
depends very slightly on ( in the region < ( < 0.25 (Fig. 18). In this case the clusters 
with diameter ~ 7ia + /2 « 1.5a + are formed. Each cluster is composed of a large positive 
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(negative) ion surrounded by a narrow layer of compensating negative (positive) charges. 
Such clusters are associated with pre-transitional ordering - they are present in the solution 
when the transition to the crystal is approached. Stability analysis alone is not sufficient 
for determination of the crystalline structure. We can expect formation of the colloidal 
crystal with periodic distribution of particles surrounded by a cloud of counterions for very 
small volume fractions. It should be noted that at the point of the A-line corresponding 
to C = 0.05 the inverse screening length is na + ~ 0.7. The screening increases with the 
increase of (, indicating the reduction of the long-range repulsion along the A-line. For 
0.1 < £ < 0.2 we find 9 > 80, and the OP consists of almost pure number density wave 
with the wavelength fa 1.5cr + , with negligible periodic ordering of the charge density (see 
Fig. 19). We may thus conclude that nearly charge-neutral units composed of the particle 
surrounded by a thin layer of the neutralizing counterions are formed. Crystal formation of 
neutral (thus noninteracting) units is not expected if their volume fraction is ( ~ 0.1. Thus, 
a re-entrant melting could be expected in this region. For ( > 0.25 we find decreasing 9 and 
the OP consists of both the number and the charge waves. The period of the OP oscillations 
increases with the increase of the volume fraction. Such a behaviour may indicate the ion 
rearrangement leading to the formation of larger clusters. 

In experiments a fluid-bcc crystal-fluid phase coexistence was found with an increase 
of the colloid volume fraction Q|- The dilute fluid - dilute crystal - dense fluid - 
dense crystal transitions, for increasing £, with the re-entrant melting for £ ~ 0.1 were also 



seen [ll| . In the latter experiment room temperature in our reduced units is T* fa 0.39. For 
0.1 < Q < 0.2 the reduced temperature at the A-line is much higher than 0.39, and further 
studies are necessary for comparison between our theory and experiment. Nevertheless, 
we can make an observation that crystal phases were observed for the volume fractions 
corresponding to the OP consisting of both the number and the charge density waves, and 
the re-entrant melting was observed when the OP consists only of the number density wave, 
with no charge inhomogeneity. 

V. CONCLUDING REMARKS 

In this paper we have used the CVs based theory to study the periodic ordering in the PMs 
with size and charge asymmetry. We consider the Gaussian approximation of the functional 
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of the grand partition function which, in turn, leads to the free energy and the direct 
correlation functions in the RPA. Using analytic expressions for direct correlation functions 
in the RPA we study the effects of the size and charge asymmetry on the instabilities of the 
uniform phase with respect to the periodic ordering. 

We determine the CV associated with the OP. To this end, we diagonalize the square 
form of the Hamiltonian and analyze the behaviour of the eigenvalues. This enables us to 
identify the OP connected with the periodic ordering and determine the character of the 
dominant fluctuations along the A-lines. 



We extend the study undertaken in Ref. 



by introducing several modifications. We 



take into account the dependence of the reference hard-sphere correlation functions on the 
wave vectors and adopt the WCA regular izat ion of the Coulomb potential inside the hard 
core. While the latter leads to a significant decrease of the A-line temperatures, the former 
gives rise to their increase for the same regularization scheme. The combination of the 
above-mentioned modifications gives rise to the A-line temperatures much lower than those 
found in Ref. j^sj. A similar situation holds for the whole range of diameter-, A, and charge, 
Z, ratios of the two ionic species. Besides, the WCA regularization scheme leads to the 



higher values of t 
compared to Ref. 



le wave vectors kf, (smaller = 2n/k} J ) associated with the A-line when 
. Precise values of kb cannot be determined in a perturbative approach, 



since the regularization of the Coulomb potential is not unique. 

According to a typical behaviour of the A-lines, we distinguish three regimes: the regime 
of small size asymmetry regime (A < 2), the regime of moderate and large charge asymmetry 
(2 < A < 20) and the regime of very large size and charge asymmetry (A = v — 200). As 
is seen, the first regime is found to be narrower than that considered in the previous work. 
Remarkably, the qualitative dependence on the charge asymmetry appears already in the 
first regime: for a small charge asymmetry, the A-lines T* versus ( demonstrate the departure 
from the monotonously increasing behaviour which becomes more evident with an increase 
of A. This change in the phase diagram is directly related to the non-local approximation 

our results also show that the 



used for the reference hard-sphere system. Unlike Ref. [28 1 
models with A > 4 and Z < 1 (a large charge at the smaller ion) do not undergo the 
instabilities with respect to the periodic ordering in the whole region of (. Thus, we can 
state that the modifications introduced in this paper have led to quantitative and partly 
qualitative changes in the phase diagram obtained in Ref. 
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We conclude that both the size and charge asymmetry affect the periodic ordering in the 
systems with dominant Coulomb interactions, namely: (1) for A < 2, the charge density 
oscillations dominate in the OP - the phase transition to an ionic crystal with a compact 
unit cell is expected; (2) for A > 2, both the charge density and the total number density 
oscillate along the structural line and their contributions to the OP depend on A and Z\ 
(3) for a moderate and large A, a crystalline phase can be stable for intermediate volume 
fractions and re-entrant melting occurs at high volume fractions, the large nearly neutral 
clusters may be formed at small volume fractions; (4) for a large A and small Z, the fluid- 
crystal phase transition can be preempted by the gas-liquid-like phase separations; (5) for a 
very large A and Z, the colloidal crystal with periodic distribution of particles surrounded by 
a cloud of counterions is expected for very small volume fractions and a re-entrant melting 
for the higher volume fractions. 

The RPA predicts only the existence of a region in which a model ionic fluid is unstable 
with respect to the OP oscillations associated with the periodic ordering. In this context 
our phase diagrams indicate the pre-transitional effects. The fluctuation effects of the higher 
order than the second order should be taken into account in order to get the information on 
both the more precise location of the phase diagrams and the pattern shapes which can be 
formed. This will be done elsewhere. 
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VI. APPENDIX 



A. Explicit expressions for the coefficients tu 



where 



, a 22 (k) a 12 (k) 

Inn — — ^ — , tNQ — ^ — , 

a 2 i(fc) an(fc) 

l QN = a — > tQQ - — 7 — ' y Z< > 



V 1 + a\ VI + o| 

a 2 i(A;) = -7==, M*) = -7==f. ( 28 ) 



A = ana 2 2 - ai 2 a 2 i = ^^y^' , (29) 

CNQ{k) 

c QQ (k) - c NN (k) ± J(c NN (k)-c QQ (k)y + 4c NQ (k)* 

aiAh) = 2^ • (30) 
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